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1.  Introduction 

A  crystalline  solid  is  generally  considered  to  be  an  assembly  of  almost  periodically- 
spaced  atoms  or  molecules.  A  set  of  periodic  atomic  equilibrium  positions  is  postu¬ 
lated  to  exist  such  that  no  net  force  acts  on  any  atom  when  all  are  at  their  equilib¬ 
rium  positions.  Thermal  excitation  (or  quantum  mechanical  zero  point  energy  at 
low  temperatures)  causes  the  atoms  to  vibrate  about  these  equilibrium  positions. 

The  atoms  interact  with  each  other  through  the  much-studied  interatomic  forces 
so  that,  if  a  single  atom  is  displaced  from  its  equilibrium  position,  a  force  acts  on  the 
others  (and  indeed  a  restoring  force  acts  on  the  displaced  atom).  The  magnitude  of 
this  force  decreases  rapidly  with  an  increase  in  the  interatomic  distance.  It  is  gener¬ 
ally  assumed  that  the  displacements  of  atoms  from  their  equilibrium  positions  are 
so  small  that  the  interatomic  forces  between  a  pair  of  atoms  is  proportional  to  the 
deviation  of  their  separation  distance  from  its  equilibrium  value.  The  force  constant 
is  largest  for  nearest  neighbor  pairs  in  the  crystalline  lattice,  and  becomes  very 
small  for  more  distant  neighbors.  In  view  of  this  Hooke’s  law  approximation  to 
interatomic  forces,  a  crystal  can  be  visualized  as  a  periodic  array  of  coupled  springs 
and  masses.  Such  a  system  of  coupled  oscillators  has  a  set  of  normal  modes  of  vi¬ 
bration  in  terms  of  which  all  motions  of  the  system  can  be  described. 

The  theory  of  lattice  vibrations  has  become  of  considerable  importance  for  several 
reasons.  The  thermodynamic  properties  of  a  crystal  depend  on  the  manner  in  which 
its  lattice  vibrations  are  excited  by  its  thermal  energy.  The  optical  properties  of 
ionic  crystals  are  determined  by  the  character  of  the  excitation  of  lattice  vibrations 
by  electromagnetic  waves.  The  electrical  properties  of  superconductors  and  semi¬ 
conductors  seem  to  be  influenced  by  lattice  vibrations. 

In  the  past  few  years  considerable  progress  has  been  made  in  the  detailed  under¬ 
standing  of  lattice  vibrations  through  the  investigation  of  simplified  models  and 
general  topological  theorems  [  1  ]  through  [6] .  We  shall  concern  ourselves  here  with 
the  analysis  of  a  model  of  a  simple  cubic  lattice  with  interactions  between  nearest 
neighbors  only;  “cubes”  of  1,  2,  3  and  a  very  large  number  of  dimensions  will  be 
discussed.  It  is  well  known  that  noncentral  forces  must  be  included  in  this  model  if 
it  is  to  be  stable  against  shear.  The  distribution  function  of  the  vibrational  frequen- 
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cies  of  normal  modes  of  oscillation  of  this  model  has  been  investigated  by  Rosen- 
stock  and  Newell  [6]  who  revived  interest  in  the  model.  This  model  has  the  unde¬ 
sirable  property  that  displacements  of  atoms  in  the  directions  of  the  various  coordi¬ 
nate  axes  are  independent  of  each  other.  It  has  the  advantage  that  most  of  its 
interesting  properties  can  be  described  in  relatively  simple  analytical  forms,  a 
feature  that  is  difficult  to  duplicate  in  more  complicated  models. 

The  first  problem  to  be  considered  here  is  the  determination  of  the  number  of 
normal  modes  in  a  given  frequency  interval.  This  quantity  is  required  in  the  calcu¬ 
lation  of  the  thermodynamic  properties  of  a  crystal. 

The  second  problem  discussed  is  the  distribution  function  of  the  location  of  a 
given  atom  with  respect  to  its  equilibrium  position.  This  problem  has  been  examined 
by  Peierls  [7] ,  and  briefly  by  Wigner  [8] .  We  shall  show  that  this  distribution  func¬ 
tion  is  Gaussian,  and  find  analytical  expressions  for  the  dispersion  in  terms  of 
dimensionality,  temperature,  and  interatomic  forces.  In  the  early  days  of  X-ray 
crystallography,  Debye  [9]  investigated  the  effect  of  this  Gaussian  distribution  on 
the  broadening  of  X-ray  spots  or  lines. 

The  third  problem  to  be  mentioned  is  the  effect  of  local  disturbances,  such  as 
impurities,  holes,  etc.,  on  lattice  vibrations.  A  method  will  be  outlined  for  handling 
these  situations. 

Finally  we  shall  give  a  brief  calculation  of  the  quantum  mechanical  zero  point 
energy  of  our  lattice  model. 

Many  of  the  mathematical  problems  discussed  here  also  appear  in  the  theory  of 
random  walks  on  discrete  lattices  and  in  the  theory  of  the  tight  binding  approxi¬ 
mation  of  electrons  in  solids. 

Since,  in  our  model  the  x,  y,  and  z  vibrations  are  independent,  we  can  obtain  our 
required  results  by  considering  the  vibrations  of  lattices  with  one  degree  of  freedom 
associated  with  each  lattice  point. 

2.  Normal  modes,  Slater’s  sum  and  thermodynamic  quantities 

Let  us  consider  a  set  of  identical  particles  of  unit  mass  (the  mass  can  be  inserted 
in  final  formula  to  make  the  units  come  out  correctly),  each  having  one  degree  of 
freedom  and  each  being  coupled  to  its  nearest  neighbors  on  a  n-dimensional  simple 
“cubic”  lattice  with  ( N  +  2)  particles  along  each  edge  of  the  n  dimensional  cube. 
We  choose  all  the  force  constants  in  a  direction  parallel  to  a  given  “cube”  axis  to 
have  the  same  value  and  postulate  the  potential  energy  of  interaction  to  be 

N 

(2.1)  =  (w»i|  ,mj  ,mg ,  •  •  •  1  .ntj ,  •  •  •)* 

wi  j  ,m2 1  ■  • ’“O 

N 

I  (Urny  ,m2  ,*nj  >•  •  •  Wm  T"  *  *  *  • 

m[,m  2>’*’— 0 

The  constant  yi  will  generally  be  used  to  represent  the  central  force  constant,  and 
the  other  y/s  the  noncentral  ones  which  are  generally  smaller  in  value.  We  have 
chosen  to  be  the  deviation  of  the  configuration  of  the  particle  at  the 

(mi,  m2,  •  •  *,  m„)-th  lattice  point  from  its  equilibrium  value.  We  shall  choose  the 
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boundary  conditions  to  be  such  that  all  particles  on  the  cube  faces  are  fixed: 

(2.2a)  Uo  ,m3  ,mg  1  •  •  •  =  1  .mj  ,mg  ,  •  •  •  0  , 


(2.2b) 


,0  ,013,  •  • » 


,N  + 1  ,mg  ,  •  •  •  0  ;  etC. 


All  u’s  with  a  subscript  0  or  N  +  1  are  chosen  to  be  zero.  We  shall  be  interested  in 
systems  in  which  N  is  very  large,  say,  0(1023),  and  in  the  limit  as  N  -*  <» . 

Since  the  kinetic  energy  of  our  system  of  particles  is 


(2.3) 


T  = 


2 

m1  ,m2  ,m3  ,••• 


> 


the  Hamiltonian  H  =  T  +  V  and  the  Lagrangian  L  =  T  —  V  are  quadratic  forms. 
They  can  be  diagonalized  through  the  introduction  of  the  normal  coordinates 
i which  are  defined  so  that 

(2.4) 


__  /  2  V/2  f 

VV  +  1/ 


1«  »•  •  •  sm 


TlTllSl  .  XW12S2 


sm 


jy  +  i  —  i\r  +  i 

If  we  write 

N 

(2.5)  =  )  ]  (Wmj  ,»i2 .  •  ■  •  ,m.j  ,•  •  •  ) 


N  ,  N 

Tins  .  Tins  y 

_o  AT  +  1  Sm  AT  +  1  —  AT  +  1  AT  +  1  2 


and  use  the  fact  that 

TV 

(2.6)  X sin 

and 

(2.7) 

we  see  that 

(2.8)  £ 


Tins  this'  1  />T  .  . 

cos  \r" r  cos  „  ,  ■■  «■  =  -  (A  +  1)  8,,'  , 


Trms  .  Tins 


S,  C°S  iV  +  1  AT  +  1 


sin 


=  0, 


and  hence  that 

(2.9) 

where 


*  - 2  JL.0  -  ~  f+i)  <•■••• 

i  £ 

-5  E 


4>  =  jr  X,  012  7?  2  , 


2§7'(1  COSJV+l)' 


(2.10) 

Frequently  we  let 

(2.11)  <f>i  =  TSj/ (N  +  1)  ,  «y  =  1,  2,  •  .  N  . 

The  largest  value  of  «2  is 


c CL 


—  4(71  +  72  +  *  ’  •  +  7n)  • 


(2.12) 
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The  Hamiltonian  of  our  system  is,  in  normal  coordinates 


(2.13) 


H  =  k  £  re . +«£ . < . ]• 


s1s2...  =  1 


This  Hamiltonian  leads  to  the  Schroedinger  equation 


(2.14) 


£  {' 

'1*2'  ‘  ’■“1  '• 


-  ,  +  (2 Eg  a  -  co2  g  V] 

dlj"  ...  *  1  *2*  *  1  *2*  ‘  *  * 


..... .)»}  -  0 


(2.15) 


E  =  E  . 
*1*2‘  *  ■”1 


The  variables  separate.  A  typical  wave  function  with  its  associated  energy  level  is 


(2.16a) 


(2.16b) 


N  1 

■®,|”*1«2'"}  =  s  „^._1  ^  0,*i*2'"(n*l*2'"  2^ 


(2.16c) 


2  =0)  7J2  /^  . 

*1  *2'  '  '  *1*2*''  '*1*2"  *' 


The  brackets  {nSl.,...}  and  }  represent  the  sets  of  all  n’s  and  rj’s  which  refer 

to  all  states  and  coordinates  of  particles.  The  functioning)  is  defined  as 


(2.17) 


i»(x)  = 


e  1/?J  '2  Hn(x) 
(2nn  !7t1/2)1/2 


Hn(x )  being  the  nth  Hermite  polynomial. 

The  position  distribution  function  of  a  system  of  particles  at  equilibrium  at 
temperature  T  is  proportional  to  the  Slater  sum 


(2.18) 


S(u )  =  E  l*n(«)r  exp  (- EJkT )  , 


where  {i„(n)}  is  the  set  of  wave  functions  of  our  lattice  and  \En)  the  set  of  asso¬ 
ciated  energy  levels.  The  Slater  sum  of  our  system  of  interest  is  the  product  of  Nn 
factors  of  which  the  (sx,  s2,  •  •  *)-th  is 


(2.19)  Ssi,sv...(xw..)  =  (a:n«2...)exp  {-(n  +  \)hunS2.../kT} 


2t  sinh 


'1*2’"  -!/2  ( 

exp^— oigj* 


.h  tanhl- 


(l  ^*1*2"^ 

\2  kT  // 


(this  summation  is  given  by  Titchmarsh  [10]).  This  expression  was  derived  in  an 
interesting  manner  by  Bloch  [11]. 
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The  logarithm  of  the  partition  function  of  our  system  of  particles,  being  the 
logarithm  of  the  integral  of  the  Slater  sum,  is 

AT  4  1  *  ■ 

(2.20)  log  Z  =  ^2  i  log  2  \ csch  2  kT 

The  various  thermodynamic  properties  of  our  system  can  be  derived  from  log  Z. 
When  the  number  of  degrees  of  freedom  becomes  large  the  frequencies  become  dense 
and  one  can  introduce  a  frequency  distribution  function,  or  frequency  spectrum, 
g(v)  (where  2irv  =  u)  with  the  property  that  f$\  g(v)dv  is  the  number  of  frequencies 
between  vi  and  v2.  The  log  Z  is  expressed  in  terms  of  g(v)  as 

(2.21)  log  Z  =  Jo  g(v)  log  {|  sech  (|  |j;)}  dv 
where  vl  is  the  largest  frequency. 

The  first  statistical  problem  associated  with  lattice  vibrations  that  we  shall  con¬ 
sider  here  is  the  determination  of  the  distribution  of  frequencies,  or  frequency 
spectrum  of  our  model. 

3.  Frequency  spectrum 

We  have  found  the  circular  frequencies  of  our  model  to  be  given  by1 

(3.1)  wij.,...  =  2  2  7/  (l  “  cos  aT+'i)  =  4  ZTysin2 

with  Sj  —  1,2 ;  •  •  • ,  N  and  j  =  1,  2,  •  •  • ,  N. 

The  values  of  w2  are  finite  in  number,  but  as  N  becomes  large,  the  random  variable 
to2  has  a  limiting  density  Gn.  That  is,  if  we  let  (a)  be  the  number  of  choices  of  s, 
such  that  w2  ...  <  a,  then 

/a 

Gn(x)dX  . 

-co 

The  density  Gn  can  be  computed  by  the  method  of  characteristic  functions.  Con¬ 
sider  (1  /Nn)Hn  as  a  cumulative  distribution  function.  Then  the  characteristic 
function 

(3.3)  /.%)  =  E(eim) 


_1_ 

Nn 


N 


z 

fo  •  •  • 


exp  {2 ia  ^  jj(l  —  cos  tSj/[N  +  1])}  . 


In  the  limit  as  N -*■«>,  this  sum  reduces  to  the  integral 

(3.4)  /„(«)  =  7r_n  J*  •  •  •  J* exp  ^2ia  2  7y(l  -  cos  0y)|  d4>\  •  •  •  d<f>n  . 


We  have  considered  the  set  of  all  frequencies  generated  by  (3.1)  to  represent  a 

1  The  energy  levels  of  electrons  in  simple  cubic  lattices  also  satisfy  this  formula  when  one  uses 
the  tight  binding  approximation  and  considers  only  interactions  between  nearest  neighbors.  Hence 
Gn{E)  represents  the  energy  distribution  function. 
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population  in  which  each  set  of  s/s  has  the  same  probability  of  occurring.  The 
function  Gn(o>2)  is  the  probability  density  function  of  a  “large  number”  of  squares 
of  circular  frequencies  chosen  at  random  from  our  population. 

The  function  Gn(o>2)  is  related  to  the  g(v)  defined  above  (2.21)  by 

(3.5)  g(v)  =  4TQ}NnG(ui2)  . 

Since  the  Bessel  function  J0(x)  has  the  integral  representation 


(3.6) 


«T' 

3 

o 

CM  _l 

3 


Jo(x)  =  7T"1  f  e~ixcoa*d<t>  , 
Jq 


Frequency  spectrum  of  a  1  -D  lattice 


we  can  rewrite  (3.4)  as 

n 

(3.7)  Mo)  =  n  exP  (2ta7J)Jo(2a7i)  . 

i-i 

Hence  (see  also  [  12] ) 

(3.8)  <?„( u>!)  41"  fl  lM2c,y,^}da  . 

Z-rr  J  _ro  /=  1 

The  largest  frequency  corresponds  to  =  fa  =  •  •  •  =  x  and  has  the  value 

(3.9)  <1>L  =  4(71  +  72  +  •  •  •  +  7n)  , 
where  we  define  by 

(3.10)  n/3y  =  7l  +  72  +  *  •  •  +  7n  • 

We  see  that 

(3.11)  G(u)  =  0  if  w2  >  oi\  =  4n/3  ,  or  «  <  0  . 

We  shall  now  find  explicit  analytical  expressions  for  Cr„(o>2)  when  n  =  1,2,  and  3, 
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and  an  asymptotic  form  for  large  n.  We  abbreviate  the  phrase  n-dimensional  by 
n  —  D  and  the  words  frequency  spectrum  by  FS. 

(a)  n=l.  The  1  —  D  expression  for  (3.5)  is 


(3.12)  G  i(co2)  =  (2 


i,Tr'  JL 


exp  [  —  ia(o)2  —  2yi)]Jo(<2cxyi)da 


l/xco(coi  —  co2)1/2 


•  r  2-4 
II  CO  <.0)1, 


•f  2  .  2 

II  CO  >  COZ,  . 


Equation  (3.5)  implies  the  following  (see  figure  1) 

(3.13)  ^(jO  =  2W7r_1(l  —  f2)~1/2  with  /  =  v/vl  , 
where,  if  we  admit  a  nonunit  mass  M, 

(3.14)  vL  =  (2t) ~ 1  (47 1 / M) 1 7 2  . 

At  low  frequencies 

(3.15)  ff(v)  ~  2N(M/yi)v*  . 

(b)  n  =  2.  The  2-Z)  frequency  density  function  is 

(3.16)  ft(«!)  =  ~  f  J0(2a7,)Jo(2«Y2)d«  • 

When 

(3.17)  «/ 0(2x71)  J 0(22:72)  =  —  I  / o(2x[7i  +  72  —  27172  cos  0]ly 

7T  «/0 

is  substituted  into  (3.16)  and  the  formula 

/•„  f  (o,  -  6!)-,,!  if  a  >  5 


(3.18) 


«/0(aZ)  cos  &Z  dZ  = 


if  a  <  b 


is  applied,  it  is  found  that 


(3.19)  G2(o)2)  = 


°  [w2(wl  —  co2)  —  I67 172  cos2  |0]1/2 

if  o)2(o)l  —  co2)  >  I67 172 

p _ d0 _ 

fl°  [co2(col  —  CO2)  —  I67172  COS2  ^0]1/2 

if  co2(co2  —  CO2)  <  1 67 172 


where 

(3.20) 


cos2  %do  =  co2(col  —  co2)/167i72  . 


2i6 
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(r2(co2)  is  immediately  expressible  as  a  complete  elliptic  integral  of  the  first  kind 
in  both  ranges.  We  have 

(3.21)  (?2(co2)  =  - - - - K  (  4(7l7z)  ^  ^ 

coir2(col  —  co2),/2  \co(ajL  —  a)2)1/2/ 


if  oj2(col  —  co2)  >  167172  , 


(3.22)  G2(co2)  =  - - - K 

27T2(yi72)1/2  4(7i72)1/2 


if  0  <  co2(co2  —  co2)  <  167i72  . 


0  .2  .4  .6  .8  1.0  0  .2  .4  .6  .8  LO 

*  2 

f  =  (co/cuL)  t  =  v/vL 

Figure  2 

Frequency  spectrum  of  a  2-E  lattice  with  72/71  =  1/9.  Logarithmic  singularities 
occur  at  /  =  0.316  and  0.948. 

The  inequality  associated  with  (3.21)  is  equivalent  to  (co2  —  471) (o>2  —  472)  ^  0. 
Hence  (3.21)  is  valid  when  472  <  co2  <  471  (we  assume  72  <  71).  Similarly  (3.22)  is 
valid  when  co2  >  471  or  co2  <  472.  It  is  to  be  noted  that  in  the  limit  as  y2  -*■  0, 
equation  (3.21)  approaches  the  one-dimensional  result  (3.12).  We  have  plotted 
<?2(co2)  and  g(y)  in  figure  2.  As  co  -►  0  (3.22)  becomes  (r2(co2)  ~  [4w(7i72)1/2]_1  so  that 

(3.23)  g(v)  ~  2ttvN2 (M 2/7i72)  1/2 . 

G3(co2)  has  a  logarithmic  singularity  at  co2  =  472  and  one  at  co2  =  471. 

(c)  n  —  3.  The  3-D  frequency  density  function  is 

(3.24)  <?,M  =  L  J  "  J^2eerOM2aY^M2ce,,-)da  . 

The  function  (j3( co2)  was  first  computed  by  Bowers  and  Rosenstock  [2]  for  the  case 
7i  =  72  =  73.  Other  cases  have  been  given  by  Rosenstock  and  Newell  [6].  The 
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integral  can  be  expressed  in  terms  of  generalized  hypergeometric  functions  of  three 
variables,  which  in  the  range  of  convergence  has  the  series  representation 

(3.25)  F3(a\  0i,  02,  03  j  7i,  72,  73;  21,  22,  23) 


=  E 

171  1^2  *  *  * 


^m1+m2+ms^l)m1^2)m2^3)m3 
-  7m  1  ^m2  ^m3 

(yiLSyz)rnSyz)mam'lm*]mil  1  2  3 


Here  (a)K  =  r(a  +  v)/T(a).  We  use  the  formula  (see  [13]) 


(3.26) 


f 

J  0 


e  pa  Jo(oeyi)  J0(cey 2)  Jo(ay3)da 
1 


1  F-  ( i.A  1  1.1  1  1.  ^inn 

"3  I  1  >  2 ,  2  J  2  ,  A ,  A,  1  ,  ,  — 

p  +  §«o>2  V  p  +  £iw!  p 


2t72 


+  hiuL  V  +  %io>L, 


2278  \ 

+  W/ 


to  obtain 


(3.27)  Gs(o?2)  =  ^  F,(l;  h  h  h;  1,  1,  1;  27^“2,  272a,-2,  2y*T') 


_j_  - Fs(l;  2,  2,  h)  1,  1,  1)  27i(oj1— o?2)  *,  272(0?!— o?2)  l,  273(0?!— o?2)  *)  . 

2irt(o?! — o?2) 

Hence  G3(o?2)  is  symmetrical  with  respect  to  o?2  =  §0?!.  Since  the  required  properties 
of  these  F3  functions  have  never  been  discussed  we  shall  find  it  more  instructive 
to  analyze  our  frequency  spectrum  from  a  slightly  different  point  of  view. 

We  apply  Parseval’s  theorem  to  the  integration  of  (3.24).  It  is  to  be  recalled 
that  if 


(3.28)  fly) - —  f“  F(a)eiavda  and  g(y)  =  — —  f"  <?(«>***« , 

**  —  00  /f*  \l/2  ‘'—oo 


then 

(3.29) 

Hence,  if  we  let 

(3.30) 
and 

(3.31) 
we  find 

(3.32) 


J°°  F(a)G*  (a)da  =  J 


(2tt)*/2 

f(y)g*(y)dy . 


(7(a)  =  c/o(2o7i)  exp  ia(o?2  —  271  —  272  —  273) 
F(a)  =  J 0(2072)/ 0(20:73)  , 


G3(o?2) 


=-  r 

2t rJ.„ 


/(y)y*  (y)dy , 


where /(y)  and  0(y)  are  determined  as  follows. 
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We  have 


(3.33) 


g(y)  = 


J 0(2a7i)  exp  [ ia(y  +  o>2 


27l  _  272  —  273)]da  . 


This  integral  is  of  the  same  form  as  that  of  the  1-D  frequency  density  function  and 
has  the  value 


(3.34) 


where 

(3.35) 


(2 /V)1/2[(y  +  w2  —  a>2)(wi  —  C03  —  w2  —  t/)]“1/2 

fKl/)  =  ^  if  ca2  —  o>2  <  y  <  wi  —  w2  —  to2 

10  if  (y  a)2  —  o>l  +  co2)(2/  +  co2  —  w2)  >  0  , 

<*>1  =  2(71  +  72)  ,  o)\  —  2(71  +  73)  ,  and  to2  =  2(72  +  73)  • 


We  postulate  71  ^  72  ^  73  so  that  coi  S:  o>2  ^  oj3. 

We  find  /(y)  in  the  same  manner  that  the  frequency  spectrum  of  a  2 -D  lattice 
was  determined: 

0  if  y*  >  cos  , 


(3.36) 


f(y)  = 


\ (2/ 7T *7 27 3) 1 /2K (5  [(o>  3  -  2/V7373]172) 

if  4(72  —  7s)2  <  y2  <  £03  , 
(2A)3/V3  -  y*r1/2K( 4[7372/(£043  -  y*)]1/2) 
if  y  <  4(72  -  7s)2  . 


The  function  g(y)  is  sketched  for  a  typical  value  of  u  in  figure  3  while  f(y)  is 
sketched  in  figure  4.  Since  the  integral  G3(co2)  is  the  integral  of  the  product  of  these 


The  3-D  frequency  spectrum  is  proportional  to  the  integral  of  the  product  of  g(y)  and  f(y).  The 
functions  f(y)  and  g{y)  are  sketched  here.  f(y)  is  independent  of  «  while  g(y)  moves  to  the  left 
as  w2  increases. 


two  functions,  only  that  range  of  y  for  which  neither  f(y)  nor  g(y)  vanish  contributes 
to  G3.  Notice  that  the  nonvanishing  range  of  g(y)  is  to  the  right  of  co2  when  o>2  |  0, 
and  is  to  the  left  of  — w2  when  002  ^  o>2.  Hence  G3  =  0  when  w2  ^  0  or  ^  o^asis 
required  by  equation  (3.11). 
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Figure  5 

The  shaded  area  in  this  figure  represents  an  example  of  the  overlap  of 
f(y)  and  g(y)  when  yi  =  Y2  =  73  =  y. 


0  y  2y'SyAy  -2y-y  0  y  2y  Zy  4y 

y  y 

Figure  6 

Variations  in  f{y)g(y)  with  «2.  The  area  under  the  various  fg  curves  is  proportional 
to  the  frequency  distribution  which  corresponds  to  the  appropriate  value  of  o>*. 
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When  co2  is  very  small  a  slight  nonvanishing  overlap  of  f(y)  and  g(y)  occurs  in 
the  neighborhood  of  y  =  <x>\.  The  main  contribution  to  G3(co2)  then  comes  from  the 
1  -D  type  of  peak  in  g(y)  at  y  =  u>\  —  «2.  Hence  in  the  limit  as  w  -*■  0  we  replace  the 
elliptic  integral  (3.36)  by  its  asymptotic  value  \k  and  we  replace  the  factor 
(<j)\  —  &>2  —  co2  —  y)  in  (3.34)  by  co2L  —  2co2  =  4?i  to  find 

fo,2 

(3.37)  GzisJ1)  ~  [87t2(7i7273)1/2]_1  I  dy/(y  +  a>2  —  a>2//2 

J  u2-0>2 


=  03 / 4ir2  (7 17273) 1 7  2  • 


0  .2  .4  .6  .8  1.0  0  .2  .4  .6  .8  10 

t2  =  (w/ojl)2  f  *  &)  /wt 


Figure  7 

Frequency  spectrum  of  a  3-D  lattice  with  71  =  72  =  73  =  7.  Here  wt  =  127. 

Hence  the  frequency  distribution  function  approaches  [see  equation  (3.8)  ] 

(3.38)  ^(1^)  ~  47n/W3(M3/7i7273)1/2  as  v~^0. 

Since  Gr3(w2)  =  G3{o3\  —  «2)  in  our  model,  we  also  have 

(3.39)  g(v )  ~  Sv^ltN5 (Ms/ 7 1727 3) 1 7 2 [  1  —  (v/vLf]1/2  as  v  vL  • 

This  shows  that  g( v)  has  a  vertical  tangent  at  v  —  vl  as  has  been  predicted  by 
van  Hove  [4]. 

In  order  to  get  a  qualitative  picture  of  the  entire  FS  let  us  first  consider  the  case 
Tl  =  =  y3  =  7.  Then  the  two  peaks  in  figure  4  coincide  at  y  =  0  (see  figure  5). 

If  we  slide  g(y)  to  the  left  along  the  y  axis  (this  corresponds  to  increasing  co2)  as 
sketched  in  figure  5  the  nonvanishing  range  of  y  in  the  product  fg  increases  (see 
figures  6a  and  6b)  and  hence  C?3( co2)  increases.  Furthermore,  the  peaks  A  and  C 
come  together  so  that  G3( co2)  increases  very  rapidly.  This  rapid  rise  stops  abruptly 
at  co2  =  47  for  two  reasons.  First,  the  end  of  the  nonvanishing  part  of  g(y)  passes 
the  end  of  f(y)  at  47.  Hence  the  total  length  of  the  nonvanishing  y  range  of  f(y)g(y) 
stops  increasing  linearly  with  co2  and  remains  constant.  Secondly,  the  peaks  A  and 
C  no  longer  reinforce  each  other.  As  a  moves  further  to  the  left  (figure  6c)  the  de¬ 
crease  in  the  contribution  of  values  of  y  in  the  neighborhood  of  A  is  compensated 
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by  an  increase  from  the  neighborhood  of  B  [  B  gets  multiplied  by  the  rising  part  of 
f(y)  ]  •  On  this  basis  G3(to2)  remains  quite  constant  until  the  point  B  passes  C  (when 
a)2  =  87)  and  A  passes  E.  Then  both  the  effective  range  in  y  starts  to  decrease  and 
the  reinforcement  of  B  and  C  diminishes.  Hence  (j3(w2)  drops  rapidly  at  w2  =  87. 

The  complete  G3( co2)  curve  is  plotted  in  figure  7a,  while  that  of  g(v)  is  given  in 
figure  7b. 

The  above  argument  is  immediately  generalized  to  the  case  73  =  72  <  71.  Here 
the  distance  between  A  and  B  in  g(y),  being  471,  exceeds  that  between  C  and  D  in 
f(y),  473.  Hence  although  the  G3(<o2)  curve  flattens  abruptly  at  co2  =  w2  when  A 
and  C  coincide,  it  drops  suddenly  when  B  passes  D  (at  to2  =  47*)  only  to  rise  to  a 
new  peak  when  B  comes  in  contact  with  C  at  co2  =  co2L  —  co2.  The  second  peak  is  a 
reflection  of  the  first  about  the  line  co2  =  %co2L.  Hence  the  final  G3(co2)  curve  is  of  the 
form  given  in  figure  8a  with  the  corresponding  FS  given  in  figure  8b. 


0  .2  .4  .6  .8  1.0 


f  =  u)/u)L  =  v/vL 


Figures  8a  and  8b 

Frequency  spectrum  of  a  3-D  spectrum  with  yi  =  7*  =  7  and  71  =  87 


When  73  <  72  =  71  the  2 -D  type  factor /(y)  has  two  peaks  instead  of  one  as  in 
figure  4  and  the  reader  can  easily  verify  that  the  G3(<o2)  curve  is  of  the  form  given 
in  figure  9.  In  the  most  general  cases  73  <  72  <  7i,  singularities  occur  at  co2  =  co2, 
co2,  co2,  co|  —  co2,  co|  —  co2,  and  co^  —  co2.  However,  the  detailed  shape  of  G3( co2) 
depends  on  various  other  inequalities  which  might  exist  between  the  7’s. 

The  procedure  described  above  can  be  generalized  for  the  deduction  of  a  4 -D 
frequency  density  function  from  a  3 -D  one.  The  main  result  is  that  a  new  set  of 
corners  appear  and  that  the  angles  of  approach  to  corners  need  not  be  so  steep. 
This  increase  in  the  number  of  corners  persists  as  n  becomes  larger,  until  in  the 
limit  as  n  — ►  <»,  Gn(co2)  becomes  Gaussian  under  a  wide  set  of  conditions  on  the 
force  constants. 

The  function  (1  /Nn)H%  can  be  considered  as  the  distribution  of  the  sum  of  n 
independent  random  variables.  Hence  Gn  is  the  density  of  the  sum  of  n  independent 
random  variables  and  in  general 


Gn+ 


Gn(y)Gi(Ji 


-  y)dy  . 


(3.40) 
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Since  Gn  is  symmetric  about  =  2(yi  +  •  •  •  +  t„),  this  becomes 

(3.41)  GUi(co2)  =  f"  Gn( JL,n  -  y)G1(J  -  y)dy 

**  —CO 

=  I  Gn(y  +  +  W2  —  £ 03L,i)dy  . 

**  —CO 

(d)  n  very  large.  We  shall  now  give  a  set  of  conditions  under  which  G( co2)  be- 


»*■  <«/«/ 

Figure  9 


The  frequency  spectrum  in  a  3-D  lattice  with  yi  —  yt  —  Qy  and  73  =  71 

comes  Gaussian  in  the  limit  of  large  n.  The  deviations  of  the  square  of  the  normal 
mode  frequencies  from  their  average  value  is 

n 

(3.42)  co2  —  2(71  +  72  +  •  •  •  +  7n)  —  —2  X)  Ty  cos  6}-  . 

j- 1 

This  quantity  can  be  considered  as  a  sum  of  random  variables  x ,■  =  —2y ,•  cos  0„ 
j  —  1,  2,  •  •  •,  n;  with  first  and  second  moments 

(3.43)  ^(xy)  =  —  (27 // 7r)  I  cos  6j  ddj  =  0  , 

•J  0 

(3.44)  £(x2)  =  (47?/tt)  I  cos2  dj  ddj  =  27?  . 

•J  0 

Since  the  x/s  are  independent,  the  dispersion  of  co2  from  its  mean  value  is 
Bn  =  E[{J  -  2n/3i)2]  =  £  27'  =  2n/S2  . 

j“i 


(3.45) 
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We  can  now  find  a  set  of  conditions  under  which  Gn  is  approximately  Gaussian  for 
large  n.  By  the  Lindeberg-L6vy  form  of  the  central  limit  theorem  [14],  if  E  = 
max,-_i,—.»(2y//  y/2n&z)>  then  the  difference  between  Gn  and  the  Gaussian  distribu¬ 
tion  is  small  in  the  sense  that 


(3.46) 


Gn(u2)do>2  -  J°  4>(a)2)da)2 


<  6 E 


1/4 


where  4>  is  the  density  of  the  normal  distribution  with  the  same  mean  dispersion  as 

Gn. 

In  the  special  case  71  =  72  =  •  •  •  =  7n  =  7,  &  =  72  and  u2L  =  4n7  so  that 


(3.47) 


Gn{<J)  = 


27  [7m] 


1/2 


exp  [  —  4n(/2  - 


(/  =  c o/wi  =  v/vL)  , 


which  means  that  the  frequency  density  function  approaches  a  5-function  with  its 
peak  at/2  =  \  as  n  -*■  <» .  The  frequency  spectrum  becomes 

(3.48)  vLg(v)  —  4(n/7r)1/2/iVn  exp  [  — 4n(/2  —  |)2]  . 

The  behavior  of  Gn(w2)  for  very  small  values  of  <*>  can  be  determined  by  noting 
that  small  values  of  are  associated  with  small  values  of  fa,  fa,  •  •  •  so  that 

~  W/{N  +  1)  ] 2  ^2  7 ,-Sj.  The  fraction  of  frequencies  with  w#l>...  <  w  («  small) 
is  then  proportional  to  2~n  times  the  volume  of  the  ellipsoid  ^  ( 8i/ai )2  =  «2>  with 
«/-  l(N  +  D/xb,-"2  or 

(1 v  +  lY  rrrnf  (^)in 

V  2i  /  r(l  +  J»)  (7.7!  •  •  •  7n) ' 72  ' 


The  proportionality  constant  is  (t/W  +  l)2,  the  volume  of  one  unit  cell  in  the  lattice 
whose  lattice  points  are  {tts,/N  +  1 } .  Hence  the  frequency  density  function  becomes 


(3.49) 


Gn(<o2) 


n(ir-u,na') 


2r(l  +  ^tt)(7i72 


yn)11 


as  to 


0  . 


The  reader  can  easily  verify  that  this  checks  with  the  special  case  n  =  1,  2,  3 
examined  above.  Since  Gn( co2)  =  Gn(u\  —  w2)  we  also  have 


(3.50) 


&t-1/2)Vl  -  a,2)*"-1 

r(l  -1-  in)(7i72  •  *  *  7n)1/2 


2  2 

as  to  — *  . 


We  close  this  section  with  a  few  remarks  about  the  distribution  of  characteristic 
values  of  the  n-dimensional  Laplace  difference  operator  D 2.  When  n  =  1,  Z>2  is 
defined  by  D2um  =  um+ 1  —  2wm  -f-  wm_i.  The  characteristic  values  X  which  satisfy 
D2u  +  Xw  =  0  under  the  boundary  conditions  (2.2)  are  given  by  (2.10)  with 
7i  =  72  =  •  ■  •  =1: 

(3.51)  X,lf,2....  =  2  2  (l  -  cos^-j)  ,  8j  =  1,  2,  •  •  •,  N  , 


whereas  those  of  the  corresponding  differential  operator  V2  with  V2w  +  Xw  =  0  are 
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(3.52)  XM.2...  =  +  S2  +  •  •  •  +  si),  Sj  —  1,  2,  3,  •  •  *, 

if  u  vanishes  on  the  boundary  of  an  n  dimensional  cube  with  sides  of  length  L. 

It  is  well  known  that  the  number  of  characteristic  values  X,...  between  X  and 
X  +  d\  is  X-1+n/2dX.  In  our  discrete  problem,  in  the  limit  as  N  ->  »  this  number  is 
NnGn(\ )  with  Gn(\)  being  given  by  (3.8)  and  with  7x  =  72  =  73  —  1.  As  has  been 
discussed  above  various  peaks  and  singularities  occur. 

The  manner  in  which  the  distribution  function  of  the  discrete  lattice  degenerates 
into  that  of  the  continuum  is  clear  if  we  let  N  =  aL  in  (3.51),  a  being  the  lattice 
spacing,  and  rewrite  Du  +  \u  =  0  as  ar2D2u  +  \u  =  0.  Then  (3.51)  becomes 


The  largest  value  of  X,  2/a2  approaches  infinity  as  a  -*■  0.  The  characteristic  value 
which  corresponds  to  a  fixed  set  of  s/s  approaches  (3.52)  in  this  limit.  All  the  peaks 
and  singularities  recede  to  infinity  as  a  -*■  »  so  that  the  density  function  has  its 
continuum  form  proportional  to  X-1+n/2  in  any  preassigned  finite  region.  Of  course, 
in  the  limit  as  a  -*■  0,  ar2D 2  -*■  V2. 

4.  Localizability  of  particles  on  a  lattice 

We  shall  now  determine  the  distribution  function  of  a  given  particle  about  its 
equilibrium  position  in  our  lattice.  If  the  dispersion  becomes  large  compared  with  a 
lattice  spacing,  we  can  no  longer  associate  a  particle  with  a  given  lattice  point  and 
therefore  cannot  consider  a  periodic  lattice  to  exist.  This  question  of  dispersion  due 
to  heat  energy  in  a  crystal  was  first  considered  by  Debye  [9]  who  asked  if  the 
general  character  of  an  X-ray  diffraction  pattern  of  a  crystal  was  affected  by 
lattice  vibrations. 

We  can  express  the  characteristic  function  of  the  displacement  of  the 

(mi,  m2,  •  •  -)-th  lattice  point  from  equilibrium  as  the  following  integral  over  the 
Slater  sum: 


(4.1)  ..(<*)  = 


f.  r.  n 

•/  -OO  •/  «1>2'****1 

/■'■f  XI  $«1«2---(W*1»2"")^W*1*2 

— 00  **  *1  *2  *  *  * 


In  view  of  (2.4),  (2.19),  and  the  formula 

(4.2)  P°  eibve-av*dv  =  fr/a)*  exp  (~b2/4a)  , 

— 00 

we  find  that  wmim*...  has  a  Gaussian  characteristic  function 


(4.3) 


•  ’{pi)  —  exp  (jlTnipj. .  .0  )  , 
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and  hence  the  Gaussian  probability  density  function 

(4.4)  =  (27rO-iim2...)-i  exp  (— W2/2o-mim2...) 

where 


(4.5) 


llm2 


.  „  irSiWi  .  ,  rS2in2  .  ,  irsnmn 

n  sin2  77——  sin2  -  -  •  •  •  sin2  vr . ,  - 7 

...  =  un  Z  -JLtl _ jLLL _ !L±± 

hut , ... 


*1*2'  *  ‘  sn=1 


a}glS2“.tanh  2 


kT 


Since  we  are  concerned  only  with  limit  results  that  are  appropriate  for  very  large 
lattices  we  let  N  -►  <»  and  find  the  following  integral  representation  for  <r^im2... : 


(4.6) 


=  Lf.:.r 

W  wnJ  „  J 


XI  sin2  <l>jmj  <20i  •  •  •  <20„ 

/- 1 


«(0 1 


0n)\ 


The  author  has  been  unable  to  express  this  multiple  integral  as  a  simple  function. 
However,  high  and  low  temperature  expansions  are  obtainable  without  too  much 
difficulty. 

(a)  High  temperature  expansion.  When  x  is  small  (x  <7 r), 

(4.7)  coth  ^  =  -(1+3®  “  45  x  +  945  “  ‘  * 

Hence  at  high  temperatures  with  hvi/2kT  <  r, 

(4.8)  2ir[s"’1L!...+ 1  Q^)’  ~  (~)‘si,1L2+  •  •  •  ], 


where 


(4.9)  S, 


r(°)  _  -»  C  *  fsh 

mlm2’  “  X  J  ”  J  ~ 


sin2  <t>\m\  sin2  0277i2-  •  •  sin2  0„77in 


in2 


<20i<202 '  •  *  d<l>n  , 


2  2  7/(1  ~  cos  0,) 

(4.10)  Snlmr..  =  7r_n  J**  •  *  J*1 sin2  sin2  •  ■  •  sin2  0„777„<20i  •  •  •  <20n  -  (|)n, 

(4.11)  =  27r“n  J"«  •  •  Jsin2  0imi  sin2  0277i2  •  •  •  sin2  <t>nmn 

•  (  2  7/(1  ~  cos  0,)j 


<20n 


n  /*t 

=  27r“n(^x)n_1  X)  7/  I  (1  ~  cos  0)  sin2  77i0<20 
j—1  Jo 


=  (ir*  E-ri  =  (jr’-i,  etc. 

/-i 
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The  integral  can  be  reduced  to  quadratures  by  noting  that  Z~x  =  f™e~Zx  dx 

so  that  it  can  be  written  as 

(4.12)  SSm...  =  dx  n  £  sin!  *‘mi  . 

Since 

Cr  i  r 

(4.13)  I  sin2  nuf>  exp  ( xy  cos  <f>)d<l>  =  -  I  (1  —  cos  2 m<t>)  exp  ( xy  cos 


=  \  Uo(xy)  -  hm(yx)]  , 

we  have 

(4.14)  siv..  =  2’(b+i)  r*-*'#  n  -  wzr, •))*  • 

‘'O  /-I  7 

The  form  of  these  integrals  is  so  sensitive  to  the  value  of  n  that  we  shall  not  write  a 
general  expression  for  S(0)  (even  though  one  can  be  written  in  terms  of  generalized 
hypergeometric  functions).  We  shall  rather  find  S(0)  when  n  —  1,  2,  3. 

(i)  Linear  chains ,  n  =  1.  Here 

(4.15)  =  if”  <T*T  [/„(*?)  -  h„(xy)]dx 


J_ 

4y 


lim 

p-*i 


[f 


Io(x)dx 


-f 

•/ o 


e  hm(x)dx 


=  lim  (p2  -  1)  *  [1  -  (p  +  vV  -  1)  2m] 

47 

=  w/27  . 

Hence  the  one-dimensional  high  temperature  expansion  for  a2m  for  a  particle  of 
mass  M  is 

(4.16)  Gm  =  \}m  +  3  (2^)  “  90  (2^)  +•••]• 

This  series  converges  when  Kvl/21cT  <  r.  If  we  choose  m  to  be  larger  than  104  all 
temperature  dependent  terms  in  the  bracket  can  be  neglected  and 

(4.17)  al  =  if  h,L/2kT< t. 

(ii)  Square  lattice ,  n  =  2.  Here 

(4.18)  Si%...  =  |  f”  e-"r'*1'!>[/o(xTl)  -  /„l(!ry1)][J,(!ryO  -  h*,J.xy,)]dx  . 


The  formula  (see  [  13] ) 
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(4.19)  /M,„  =  f  e~px  l2»(xyi)hli(xy2)dx 

Jo 


7i'172vr(^+/i+y)r(i-f’M+j') 


Tl/2pl+2M+2 -r(1+2M)r(1+2„) 


1+H+V,  2/i+l,  2^+1,  71/p2,  72/p2) 


where  JF4  is  a  generalized  hypergeometric  function 


(4.20) 


Stia,  0;y,y';x,  y)  =  X) 


(oQm+»(fl) 


(y)m(y') nm  In ! 


m+n  in  n 

«  y 


with  (a)„  =  r(a  +  ri)/T(a)  allows  us  to  write 


(4.21) 


Sm 


lm2 


—  ft  {^0.0  Io,2m2  “f“  ^2mj,2»»2  } 


Since  asymptotic  formulas  for  generalized  hypergeometric  functions  in  the  limit 
of  large  n  and  v  have  not  been  discussed  in  the  literature,  we  shall  write  S£?mt  as  a 
form  more  appropriate  for  the  range  of  large  mi  and  m2.  We  can  determine  the 
qualitative  behavior  of  cr^imi  for  points  far  from  the  boundaries  of  our  square  lattice 
by  letting  mx  =  m2  —  m.  Then 


(4.22) 


[h(xyx)Io{xyd  -  I2m(xyi)l2m(xy2)]dx 


+ 

+ 


Jo 

I  e~x(yi*72>  Iim(xys)[Ism(%yi) 
Jo 


-  Io(xy2)]dx 

-  /0(x72)]da;] 


The  asymptotic  behavior  of  the  second  two  integrals  is  discussed  in  appendix  I 
where  it  is  shown  that  they  are  negligible  compared  with  the  first  when  m  is  large. 
The  first  integral  converges  but  if  it  is  separated  into  the  difference  of  two  integrals 
each  of  those  does  not.  However,  if  we  introduce  the  integrating  factor  exp  (— ye) 
both  converge.  Then  we  have  as  the  value  of  the  difference  (see  [  13  ] ) 


(4.23)  lim 
€*♦0 


=  lim  - 

«-*o  x 


{J\—  h(xydh(xy,)dx  -  j]” 

1  r q  ( (71+72+e)2— 71—72^  _ 
(7172) 1/2  \  27172  / 


n  {(y,+y,+tf-yl-y!)l 

2yiT2  /J 


where  the  Q’ s  are  Legendre  functions  of  the  second  kind.  Since  (see  [15])  Q,(Z)  ~ 
— §  log  (§«)  —  7  —  ¥( v  +  1),  where  Z  =  1  +  e  and  e  -*■  0,  7  =  0.57721,  and 
^(n  +  x)  —  -*(x)  =  1/x  +  l/(x  +  1)  +  •  •  •  +  l/(x  +  n  —  1),  our  integral 
becomes 

W'24)  T(7m)‘'2  I'  +  §  +  I  +  ’  ’  ’  +  2^i}  ~  '0g  2m  ■ 

Hence  as  m  -►  «> 
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(4.25)  SL°l  ~  [87r(7i72)1/2]_1  log  2 m  . 

If  we  introduce  a  mass  M  for  the  particles,  the  term  (7i72)1/2  becomes  (7172 /Af2)1/2, 
v\  =  (71  +  72) /7 r2M  and  our  desired  dispersion  becomes 

0m  2  kT  /  1  f(7i+72)2~l*  1  0  ,  1  / hvL  V  1  AvlV  ,  \ 

(4.26)  ~  2Mj,2  j27r3  [  J  log  2m  +  3  ^2A.rJ  90  \2kT)  +  *  *  * / 


where  the  series  converges  if  hvi/2kT  <  t.  Hence  as  long  as  log  2m  »  tc  (m  »  12) 
we  can  neglect  the  temperature  dependent  terms  in  the  bracket  to  find 

o  kT 

(4.27)  o-rn.m  ~  47r(7lTa)^  log  2m  hvL/2kT  <  m  . 

(iii)  Simple  cubic  lattice ,  n  =  3.  In  this  case  S^/m2tnt  is  a  sum  of  four  integrals  over 
products  of  the  Bessel  functions  [see  equation  (4.14)].  In  the  limit  of  large  mh  m2, 
m3  (lattice  points  far  from  the  boundaries  of  the  cube)  one  can  apply  the  asymptotic 
results  of  Appendix  III  (equation  III— 6)  to  find  that  those  terms  with  at  least  one 
large  subscript  become  negligibly  small.  Then  one  finds 

(4.28)  &!.%,„  2-  f"  UMUMUMdx  =  Sm 

J  a 

as  mi,  m2,  m3  ->  00 .  The  case  of  main  physical  interest  is  one  with  two  of  the  y’s 
equal,  say  72  =  73.  Then  one  finds  from  equation  11-20  and  equation  11-17  that 
Smim2m3  is  independent  of  the  m/s  as  they  become  large  and  that 

sm  =  &!«*.,  =  1(1  +  t)!  -  (7  -  D'W'W 

(4.29) 

where  7  =  (37i  +  472)/7i,  fa  —  §[  (7  —  1)1/2  —  (7  —  3)1/2]  [  (7  +  1)1/2  —  (7  —  1)1/2] , 
and  K(k3)  corresponds  as  usual  to  the  elliptic  function  of  the  second  kind  and 
K'(fa )  =  A([  1  —  /c2]1/2).  The  integral  1(a)  is  discussed  in  Appendix  II  and  plotted 
in  figure  13. 

When  noncentral  forces  are  weak  (when  71 »  72),  equation  (4.29)  simplifies  to 

(4.30)  21/2~  1)  as  yi/yt  — *  0  . 

This  result,  being  independent  of  (mi,  m2,  m3),  implies  that  at  all  lattice  points  far 
from  the  boundaries  of  our  cube 

(4.31)  °  =  {[(7i+27.)/72U(7i/7!)+  |  (hvL/2kT)'-~  (hvL/2kTf  +  --j- 

This  expression  differs  from  the  analogous  1  -D  and  2 -D  results  in  that  it  depends 
only  on  the  temperature,  particle  masses  and  force  constants. 

We  now  proceed  with  the 

(b)  Low  temperature  expansion.  When  x  is  large  coth  x  =  1  +  2e~2x  -J-  2e-4x  + 
2e_6x  +  •  •  •  so  that 
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(4.32)  vmim2...  =  ^  f  -  fa)  1  £l+2  £  exp  j  J  n  sin2  v  •  -d<£„. 

The  asymptotic  low  temperature  limit  can  be  reduced  to  quadratures  through  the 
introduction  of 


(4.33)  z~h  =  (2/tt1/2)  r  dz 

Jo 


and  a  repetition  of  the  argument  used  in  the  derivation  of  (4.14).  We  then  have 
h(i 1 


(4.34) 


2 

2* 


w*  p 

2»+i  ^0 


ar!e-*!r 


>'  n  f/.(*ry)  -  I  2m,-  (z7>)}  da; 

;= 1 


where,  as  before,  Im(x )  is  the  nth  Bessel  function  of  purely  imaginary  argument. 

We  shall  now  find  «,,•••  for  1,  2,  and  3  dimensions. 

(i)  n  =  1,  one-dimensional  lattice. 

(4.35)  cl  ~  \h  0*  f“  *-*  <T't  [/o(*y)  -  7».(»y)](fc 

*rtl/2  f  {* 03 

■  x~ie~p'I-^dx- 

=  lim  (Q-j(p)  -  Q!»-i(p)| 

where  Qn(x )  is  the  nth  Legendre  function  of  the  second  kind.  It  has  a  singularity  at 
x  —  1.  The  asymptotic  expression  in  the  neighborhood  of  x  —  1  is  (see  [15]) 


r 


r*  e~vxu 


l(x)dxi 


(4.36)  Q,(x) - §  log  (x  -  l)/2  -  7  -  *(1  +  y)  +  0(x  -  1) 


where  7  is  Euler’s  constant  0.57721...  and  where  for  integral  n 


(4.37)  *(Z  +  »)-*(Z)  =  L  +  _J_+.  . 

We  finally  obtain  the  low  temperature  limit 

(4.38)  cl  =  [*(2m  +  i)  -  *(J)]ft/(2*7*) 


.1  +  l  +  l  + 


+ 


The  sum  has  the  asymptotic  value  \  log  m  when  m  is  large.  We  can  introduce  the 
particle  mass  M  and  use  this  asymptotic  relation  to  find 


(4.39)  <j m  =  \(h/MvLTr2)  log  m  . 

(ii)  n  =  2,  square  lattice. 

We  apply  (4.34)  with  n  =  2  to  obtain 


2* 


hit  *{Foo(7  1>  72)  —  Fo,2m2(Tl>  72)— Fimi,o(7l>  72)  +  F2m1.2m2(7l,  72)} 


(4.40) 
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(4.41) 

(4.42) 


Fa  Ay i,  72) 


-r 


x~i  e-<w  uxy^ijxyi)dx 

=  (w*T,/4  r(a  +  fi  +  1)P£.  (r2)p;_#j  (r.) , 
r,  =  [(,ti  +  yi)/yj]m,  j=  1,2; 


ij  =  log2(y/y2) 

Figure  10 

Variation  of  the  integral  «r,F0o  (7172);  (see  equation  4.41)  with  (71/72)  =  2n. 


and  P'(Z)  is  the  Legendre  function  of  the  first  kind.  In  the  special  case  a  =  p  =  0, 
we  have 

(4.43)  F 00(71 »  72) 

=  tt3/2(7i72)1/4  [(ih  +  i)(r2  +  i)]1/2^  Kr*_1)(W7i5  ^t(ri-i)(7i/72)  ] , 

where  as  usual  K(k)  is  a  complete  elliptic  integral  of  the  first  kind.  We  have  plotted 
this  function  in  figure  10. 

It  is  shown  in  Appendix  IV  that  when  yXlm\  +  7 2  ’wij  is  large 

(4.44)  F 2m  j  ,2m 2(7 1 )  7*)  =  7  (  )  /(Yl^l  +  72  ^2)'  . 

4:  \7l727T/ 

Hence  as  Wi  and  m 2  , 

2  2*  for""*  „  ,  \ 

(4.45)  ~  -g  ^00(71,  72)  , 

independently  of  mi  and  m2.  Furthermore  when  72  «  71,  F0 0  has  the  asymptotic  form 

(4.46)  F0o(yi,  72)  ~  (2tt7i)"J  log  (6471/72)  . 

Hence  the  dispersion  gets  large  as  72  — ►  0. 

(iii)  3 -D  lattices. 
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The  3-D  expression  for  the  low  temperature  <7 2  which  is  analogous  to  (4.40)  con¬ 
tains  six  integrals  of  products  of  three  Bessel  functions.  Those  integrals  with  Bessel 
functions  of  order  m  1,  m2,  or  m3  can  be  expected  to  approach  zero  with  increasing 
mi,  m2,  and  m3  even  faster  than  those  of  products  of  two  Bessel  functions.  Hence 
for  those  atoms  far  from  the  boundaries,  we  have  the  asymptotic  low  temperature 
dispersion  (see  equation  (4.34)) 


(4.47) 


=  2^_  r° 

16tt*  Jo 


x~i  e-*(71+72+73 >I0(Xyi)I0(Xy2)I0(Xy3)dx  . 


A  series  expansion  can  easily  be  obtained  for  this  integral  when  72  =  73  (it  is  to  be 
recalled  that  72  «  71  in  interesting  physical  cases).  We  have 


(4.48)  4  =  ^  f"  *-»  ± 

107rV  0  m-0 


(^72)2mr(2rn  +  1) 

m![r(l  +  m)]8 


dx  . 


After  interchanging  the  order  of  integration  and  summation  and  applying  a  formula 
on  p.  196  of  Tables  of  Integral  Transforms  [13]  we  find 


(4.49)  <r20  = 

(4.50) 


_ k _ /p  (7\  1  'ST'y-Smf  _72_V  r(4m+l) 

16[72(7i+72)]1/W1/2  \  m+f  \7i+72/  [r(m+l)34 

Z  —  |(7i  +  272)/[72(7i  +  72)]1/2 


P2m-i(Z)l 


where  Pn(Z )  represents  the  nth  order  Legendre  polynomial.  When  m  is  large 
(see  [16]) 


(4.51) 
and 

(4.52) 


P  2m— 1 


(  71  +272  \ 
\2[72(7i+72)]V 


1 

(T7i)*(4m— l)1 


O’?) 


m 

[72(71+72)] 


1/4 


r(4m)/[r(m  +  l)]4 


(27r)1/228m 
8 7r2m6/2  ' 


Hence  the  mth  term  approaches 


(4.53) 


[72(71  +  72)] 1/4 
27r2m2(27i)1/2 


so  that  the  series  can  be  expected  to  converge  fairly  rapidly. 

A  good  approximation  can  be  obtained  for  <r2  by  using  (4.53)  when  m  ^  4  and 
using  the  exact  terms  when  m  <  4.  We  have 


(4.54)  X)  ™-2  =  +76)  -  (49/36) 

4 

and 

(4.55)  P_f  (2^^)  =  “[(71+72)*— 72] [72(71 +72)]  1/4A([(7i+72)2  — 72]2/7i) 

(4-56)  P,  -  l  , 

where  K  and  E  are  elliptic  integrals  of  the  first  and  second  kind.  The  Legendre 
functions  P 3/2,  Pm,  and  P11/2  can  be  expressed  as 
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(4.57) 

PS/J(Z)  =  |zP1/s(Z)-|p_,(Z) 

(4.58) 

P7/i(Z)  ■ 

=  (384Z»  -  208 Z)Pm(Z)  -  ^  (96Z!  _  25 )P.,(Z) 

(4.59) 

Pni*(Z) 

=  1(122,8802*  -  129,024Z3  +  25,668 2)P1/!(Z) 

-  (30,720 Z4  -  23,616Z2  +  2025 )P_j(Z)}  . 

We  have  plotted  the  asymptotic  value  of  alMll2<r\,2/h  as  a  function  of  72/71  in 
figure  11  as  the  temperature  approaches  0 °K. 


V  *  lo9* 

Figure  11 

Variation  of  the  dispersion  <r02  of  a  given  atom  (in  a  3 -D  lattice)  from  its  equilibrium  position 
at  absolute  zero  temperature  as  a  function  of  the  ratio  of  central  to  noncentral  force  constants 
(71/72).  The  parameter  n  is  chosen  so  that  71/72  =  2". 


In  the  limit  as  72/71  -*■  0,  the  parameter  Z  [see  (4.50)]  approaches  Kt2/yi)1/2  so 
that  (see  [15]) 

r(2B)(Ti/T,r1<< 


(4.60)  P*n-i(Z) 

so  that  the  series  in  (4.43)  becomes 


%T(2m  +  \) 


(4.61)  W TO*'4 1  2-  =  °-°666  My'r 


while  P-w(Z)  becomes  -  (  — )  log  (16  71/72).  Hence,  in  the  limit  as  71/72  -►  ® 
r  \7i/ 

and  T  -►  0 
(4.62) 


1/4 


2 

(7  0  ~ 


167t(7iM)1/2 


log  (I671/72)  . 


This  differs  from  the  corresponding  2 -D  expression  (4.46)  by  having  the  factor  16 
rather  than  64  in  the  logarithmic  term. 

A  good  approximation  to  the  first  temperature  dependent  term  in  the  expansion 
of  a2  [  see  equation  (4.32)  ]  as  a  power  series  in  ( kT )  can  be  obtained  as  follows . 
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In  the  limit  as  T  -*■  0  the  main  contribution  of  the  integrand  of 


(4.63)  t  3  JJJ &T1  23  exP  (—jhta/kT)d4>id<t>2d<t>z 

0  3 

comes  at  small  value  of  co.  However,  the  relation  between  w  and  the  <j>}s  in  this  range 
3 

is  co2  ~  ^2  y&\  Since  each  of  the  exponentials  is  a  rapidly  decreasing  Gaussian 

j=i 

function  in  this  approximation,  the  upper  limits  of  integration  can  be  extended  to  <» 
and  one  can  introduce  spherical  polar  coordinates  (after  letting  x  =  4>  17 }/2,  etc.). 
Then  the  integral  becomes 

4'64)  TWO'"  S  /.  TT  ^  <->»“/»>*» 


=  4(kT)2  yp  _2  2  (kT/hy 

T2h2(y 17*73) 1/2  3  (717273) 1 /2 ’ 

The  first  two  terms  in  a  low  temperature  expansion  of  the  dispersion  is  (in  the 
symmetrical  lattice  with  72  =  73) 


(4.65) 


2  _1_ 
00  + 


2h 


3(71  MY'2 


+ 


We  see  that  as  the  effect  of  noncentral  forces  diminishes  (that  is,  as  72/71  -*■  0)  the 
zero  point  dispersion  increases  and  the  temperature  dependent  dispersion  becomes 
effective  at  lower  temperatures.  The  next  order  terms  in  the  expansion  would  come 
from  the  4>A  terms  in  the  expansion  of  «2. 

The  dispersion  can  be  obtained  from  an  integration  over  the  frequency  distribu¬ 
tion  function  at  all  temperatures 


(4.66) 


16xW3 


f 


g(y)dv 


v  tanh  (7. 


1  hv  \ 


2  kTJ 


We  can  summarize  the  results  of  this  section  as  follows.  Every  particle  in  our 
lattice  has  a  Gaussian  distribution  about  its  equilibrium  position.  The  value  of  a2  in 
the  distribution  depends  on  temperatures,  force  constants,  and  dimensionality. 

In  a  linear  chain  <r2  is  proportional  at  high  temperatures  to  the  distance  of  the 
particle  of  interest  from  a  fixed  end  of  the  chain;  at  low  temperatures,  to  the 
logarithm  of  that  distance.  Hence  those  atoms  far  from  the  ends  of  a  long  chain 
might  vibrate  over  distances  long  compared  with  the  lattice  spacings. 

A  particle  on  a  2 -D  lattice  has  a  a2  which  is  independent  of  its  equilibrium  position 
at  low  temperatures  provided  that  it  is  far  away  from  lattice  boundaries.  At  high 
temperatures  a 2  is  proportional  to  kT  times  the  logarithm  of  its  distance  from  the 
boundaries.  Hence,  at  low  temperatures  a  particle  is  localized  at  a  given  unit  cell, 
while  in  a  system  of  a  sufficiently  large  number  of  particles  a  particle  sufficiently  far 
from  a  boundary  may  at  high  temperatures  lose  its  localization. 

The  particles  in  a  three-dimensional  lattice  with  harmonic  forces  are  localized  at 
all  temperatures.  We  have  plotted  <j 2  (figure  12)  as  a  function  of  temperature  for 
typical  particles  far  from  boundaries.  Curves  are  given  for  various  values  of  72/71. 
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As  is  to  be  expected  the  dispersion  increases  as  the  strength  of  the  noncentral  force 
constant  diminishes. 

It  is  to  be  recalled  that  in  our  nearest  neighbor  interaction  model  the  above- 
discussed  values  of  a2  correspond  to  the  components  of  displacements  from  equi¬ 
librium  in  the  direction  of  only  one  of  the  crystal  axes.  The  total  dispersion  is  the 
sum  of  that  over  all  three  directions.  When  more  distant  neighbor  interactions  are 
included  the  displacements  in  different  directions  are  correlated.  In  a  highly  aniso¬ 
tropic  lattice  the  dispersions  from  equilibrium  are  greatest  in  the  direction  of 
weakest  intermolecular  forces. 


Variation  of  the  dispersion  <r*  with  temperature  in  a  3 -D  lattice 

The  joint  distribution  function  of  the  displacements  of  a  pair  of  particles  from 
equilibrium  is  the  double  Fourier  transform  of  the  characteristic  function 

(4.67)  G({m),  {m  +  n}’,ai,  a 2)  =  E(exp  —  t[aiWmi.m2,...  +  a2Wmi+Ml,...])  . 


Here  { m }  =  (mi,  m 2,  •  •  mn)  represents  the  coordinates  of  one  of  the  particles, 
and  {m  -f-  m}  those  of  the  other.  Various  values  of  the  displacements  are  weighted 
by  the  Slater  sum.  One  finds  the  characteristic  function  (4.67)  to  be  Gaussian  with 
the  formula 

(4.68)  exp  [—  ^(ahal  +  2raia2crm<7m+M  +  al ov,+m)] 

where  the  o-’s  are  given  by  (4.6)  and 


n 

n  ( cos  —  cos  0/(2wy  +  My) } 

1 


« tanh 


(4.69)  ramcr, 


d<f>  1  d<fin  • 
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Here  r  is  the  correlation  coefficient  between  particles  separated  by  the  vector  ju. 
Since  the  characteristic  function  (4.68)  is  Gaussian,  its  Fourier  transform  is  also 
Gaussian  and  the  joint  distribution  function  is 


(4.70) 


_ 1 _ 

2Tcrm<rm+ft(l—  r2)112 


1 

2(1  -r2) 


O’m+M 


+ 


High  and  low  temperature  limits  of  r  (defined  in  (4.69))  are  obtained  in  the  same 
way  as  they  were  found  in  the  discussion  of  c.  We  shall  merely  list  these  results  for 
the  1  -D  case  and  high  temperature  asymptotic  results  for  the  3-D  lattice. 

In  the  1  -D  case  the  high  temperature  limit  of  (4.69)  is 

(4.71)  r<rm<rm+ll  =  mkT/^Mv^  . 


This  equation  combined  with  (4.16)  yields 

(4.72)  r  =  [1  +  ( n/m )]_1/2. 


Hence,  as  n  -+■  <»  for  fixed  m,  the  correlation  coefficient  vanishes.  At  low  tempera¬ 
tures  (and  finally  large  values  of  n  for  fixed  m) 


(4.73) 


ram<7r 


0 h/MvLT 2) 


+ 


1 


12  ju  -f-  1  2  n  -|-  3 


+ 


+ 


2 n  -f-  (4  m 


~  (h/2ir^Mp  £)  log  (1  -f-  2 mu  *)  ~  hm/ir2MnvL  . 

We  combine  this  with  (4.39)  to  obtain  the  asymptotic  result  as  n  -►  <»  for  fixed  m, 


(4.74)  r  ~  2(m/n)/[  (log  m)(log  *t)  ] 1/2  • 


A  direct  application  of  equation  III-2  in  Appendix  III  gives  the  3-D  high  tem¬ 
perature  result  (which  is  independent  of  m )  when  /z  and  m  are  large 


(4.75) 


TCm^m+n  — 


kT 


87T  (7 17 2 7 3) 1 ;  2s 1 Z  2M  1 


2  -1 
=  MiTi 


.  2  -1  .  2  - 

+  M272  +  M373 


Since  am  and  <rm+M  are  both  proportional  to  (kT)112,  we  find  that  r  is  proportional  to 
$-I/2.  The  value  of  the  proportionality  constant  depends  on  the  particle  mass  and 
force  constants  but  not  on  the  temperature. 

A  quantity  of  importance  in  X-ray  diffraction  is  the  distribution  function  of 
(wm+M  —  um ).  This  is  obtained  from  the  characteristic  function  (4.67)  when  ax  = 
—  a2  =  a.  Hence  (4.68)  reduces  to 

(4.76)  exp  [  ^<x  (<rm  2ra'mo'm+it  -t-  o\»+m)J  • 


Hence  our  required  distribution  is  Gaussian  with  a  dispersion 


2 

Gm  +  p  ,m 


2 

&  m 


2roma 


m+p 


+  <r 


2 

m+p  • 


In  three-dimensional  lattices  2rcrmcrm+ft  vanishes  for  large  values  of  n  and  <rm  =  <rm+M 
so  that  in  this  range  <^+M  m  ~  2<r2m.  However,  in  one-dimensional  lattices  2 r<rm<jm+» 
might  become  large.  In  this  case  we  have 
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(4.77) 

at  high  temperatures,  and 

(4.78) 
at  low. 


O'm+M.m  =  flkT /t2Mv  L 


2tt2Mv  l 


5.  Effect  of  local  disturbances  on  lattice  vibrations 

In  this  section  we  shall  outline  a  procedure  which  can  be  used  to  discuss  the  effect 
of  local  disturbances  in  a  lattice  on  its  lattice  vibrations.  By  a  local  disturbance  we 
mean  a  foreign  atom  at  a  lattice  point,  a  hole,  etc.  We  shall  restrict  our  analysis  to 
3-D  lattices  but  the  method  is  applicable  to  those  of  any  number  of  dimensions. 

The  difference  equations  from  which  we  obtained  the  characteristic  frequencies 
of  our  normal  modes  w2  in  section  2  are 

(5.1)  'Yl{Um^  +  l  ,tn  2>m  ^  2wmj  ,m2  .mg  +  Wmj  —  1  ,m2  ,m3] 

-(“  'Y’2[Wmj  .»«2+l  2Wm  j  ”1“  WTOj.m2-l.wt3] 


~1~  TsfWmj  .m2  2Wmj.TO2.TO3  ~1~  Wmj.TO2.m3— l]  "4~  -^Do2Wm  jm2m3  0  . 

Here  mh  m2,  m3  range  through  1,  2,  •  •  • ,  N.  For  convenience  we  shall  change  our 
notation  to  let  the  m’ s  range  from  —N/2  to  N/2.  When  the  boundary  conditions 

(2.2)  are  chosen  the  normal  mode  frequencies  (2.10)  result.  Now  suppose  that  some 
masses  or  force  constants  are  different  from  the  others.  Then  the  coefficients  of  the 
certain  w’s  are  different  from  those  given  above.  Indeed,  if  we  let  D  represent  the 
difference  operator  which  acts  on  wm.m,  m,  to  yield  the  above  equations,  that  is,  if 

(5.2)  DWTOJTO2TO3  =  0  , 

then  our  new  equation,  which  would  show  the  effect  of  local  disturbances,  would  be 
(the  total  displacement  of  an  atom  from  its  equilibrium  position  is  the  sum  of  the 
original  unperturbed  displacement  plus  the  solution  of  the  following) 

(5.3)  Dum jm2TO3  =  w3,k,l(mi  +  j,  m2  +  k,  m3  +  l)umi+j,m2+k,m3+i 

i.k.l 

nti,  m2,  Mz  =  0,  ±1,  ±2,  •  •  ’,  ±N / 2  , 

where  the  functions  wJ’k’l(mi  +  j,  m2  +  k,  ms  +  l )  would  characterize  the  dis¬ 
turbance.  For  example,  if  we  merely  change  the  mass  of  the  particle  at  (a,  b,  c)  to  M' 
and  leave  all  force  constants  fixed 

(5.4)  w3'k,\mi  +  j,  m2  +  k,  m3  +  l)  =  0 
unless  j,  k  and  l  are  zero  and  mi  =  a,  m2  =  b  and  m3  =  c.  Then 

(5.5)  b,  c)  =  {M  -  M')<J  . 

If  the  force  constants  between  a  single  atom  at  (a,  b,  c)  and  its  neighbors  are  changed 
all  w/s  would  be  zero  except  w;±li0,0(a  ±  1,  b,  c ) ;  w °’±1-°(a,  b  ±  1,  c) ;  w>°’0,±1(a,  b,  c  ±  1) ; 
tp°'0-°(a,  b,  c ). 
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The  set  of  equations  (5.3)  can  be  solved  through  the  use  of  Green’s  functions.  We 
shall  follow  the  ideas  used  in  the  paper  of  R.  T.  Duffin  on  discrete  potential  theory 
[17].  Since  we  are  interested  only  in  local  disturbances  let  N  °°  and  first  con¬ 
sider  the  Green’s  function  g(m  1,  m2,  ra3)  which  is  the  solution  of 

(5.6a)  Dg  =  1  for  mi  =  m2  =  m3  =  0 

(5.6b)  Dg  =  0  otherwise 

(5.6c)  g  — >  0  as  m\  +  ml  +  ml  — »  °°  . 

The  Green’s  function  can  easily  be  verified  to  be 


(5.7)  g(mi,  m2,  m3) 


1  rrr _ ei^1m1^«2^3«3)^i^2d(/)3 _ 

(2x)3  JJJ  Mu2—  2ti(1  —  cos  0i)  —  272(1  —  cos  </>2)  — 273(1  — cos  <£3) 


=  >'  da  II  j’  e 


Q~2yi  cos*i 


where  as  usual  Jm(ia)  is  a  Bessel  function  of  imaginary  argument.  This  integral 
converges  when  Mu2  >4^7,-  =  or  Mu2  <  0.  It  can  be  represented  as  a 
generalized  hypergeometric  function  of  three  variables  (see  equation  (3.27)  of  this 
article  and  [13],  p.  184).  The  function  <7(0,  0,  0)  has  been  tabulated  by  M.  Tikson 
[18]  when  71  =  72  =  73.  Since  these  integrals  occur  in  many  problems  in  which 
cubic  lattices  appear  (ferromagnetism,  electrons  in  metals,  random  walks  on  lat¬ 
tices,  etc.)  their  tabulation  would  be  a  worthwhile  project.  Asymptotic  expressions 
for  large  m’s  are  easily  obtained  by  following  the  method  discussed  in  appendix  III. 

Now  that  the  Green’s  function  is  known,  it  is  easily  verified  that 

(5.8)  umim2mz  =  s  2  ?(»ii  +  j  —  nl}  m2  +  k  —  n2,  mz  +  l  -  n3) 

jkl  mj»»2»»3 

•  w3k\ni ,  n2,  nz)unin2n 3  . 


The  application  of  the  D  operator  and  the  use  of  (5.6)  immediately  yield  (5.3).  An 
important  feature  of  this  result  is  that  it  is  exact  (no  perturbation  calculation  is 
made).  A  similar  result  has  been  obtained  by  Slater  and  Roster  [19]  in  their  appli¬ 
cation  of  Wannier  electronic  wave  functions  to  the  problem  of  an  impurity  in  a 
semiconductor.  It  has  also  been  discussed  by  Lax  [  23  ] .  The  effect  of  holes  on  lattice 
vibrations  has  been  treated  by  perturbation  theory  by  Stripp  and  Kirkwood  [20]. 

Although  several  detailed  applications  of  (5.8)  will  be  discussed  in  another  paper 
we  shall  indicate  the  approach  to  the  simple  problem  of  a  change  of  mass  of  a  single 
atom  without  the  change  in  force  constants.  This  will  give  the  spirit  of  the  method 
of  using  (5.8).  Suppose  the  impurity  atom  is  at  the  origin  so  that  w°°0(0,  0,  0)  = 
( M  —  M')u2  and  all  other  w’s  are  zero.  Then, 


(5.9) 


Wm1m2m3  =  g{mim2mz)u(M  —  M')u0oo  ■ 
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The  new  normal  mode  frequency  to2  would  be  determined  by  setting  mi,  m2,  m3  =  0 
and  solving  the  transcendental  equation 

(5.i°)  W2(M  -  M')  =  °’  0) 

for  a>2.  When  71  =  72  =  73,  Tikson's  tables  are  very  helpful  in  obtaining  the  solu¬ 
tion.  This  value  of  u>2  lies  outside  the  continuum  of  frequencies  which  we  discussed 
in  section  3.  With  this  value  of  u2,  the  vanishing  of  the  disturbance  to  the  normal 
lattice  vibrations  with  increasing  distance  from  the  impurity  can  be  discussed 
through  (5.9). 

In  cases  of  more  complicated  sources  of  disturbance  one  sometimes  has  to  solve 
a  set  of  simultaneous  equations  to  find  and  the  impurity  frequencies  a>2. 


6.  Zero  point  energy 

We  finish  this  report  with  a  brief  discussion  of  the  quantum  mechanical  zero 
point  energy  of  our  model 

(6.1)  Eo  =  %h 


(2  M) 


1/2 


h(N/rT  f-f(z  7/(1  ~  cos  *y)V,!  dpr  ■■■  dp,  . 


This  expression  can  be  reduced  to  quadratures  by  employing  the  following  represen¬ 
tation  of  x1/2: 


(6.2) 

Then 


M  2 


_J_  f°°  / 1  -  e~xat\ 

J-",  \  a2  ) 


da  . 


:(1— COS  <t> .) 


]  1 


(6.3)  Eo  -  2(27rM)1/2  C  )  /o  °  da  6 

An  integration  by  parts  finally  yields 

<«•«  -  2 %■” f  rI'2 v™  ■  ■  • 

•  [/o(6t,)  -  h(by/)]/h(by,)db  . 


d<f>n 


The  zero  point  energy  of  a  1  -D  lattice  is  [see  discussion  under  (4.34)  ] 


(6.5) 


Eo  — 


hNy 


r 

Jo 


Jj- 1/2  „-by 


e ~by  [Io(by)  -  h(by)]db 


2(2tt M)112 

-  [*(3/2)  -  ™  -  f  (zY  ■ 


Since  the  zero  point  energy  of  2 -D  and  3 -D  lattices  involves  four  and  six  integrals, 
respectively,  we  shall  not  write  them  down  explicitly  here.  We  merely  point  out  that 
each  of  these  integrals  is  of  the  form  dealt  with  in  the  discussion  of  the  asymptotic 
low  temperature  expression  of  a2  [  see  (4.34)  ] . 
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APPENDIX  I.  ON  THE  INTEGRATION  OF 


(M) 


exp  [-x(ti  +  y2)]{Imi(xyi)Im2(xy2) 


I m^ipcy  ^Im^ixy^) }  dx 


WHEN 


j  2  2  -1/2  1  2  -1/2  j  7  2 

k  1  =  WJ171  +  WI272  and  k  2 


2  -1/2  1  _ 2  -1/2 

m3yi  +  men 


ARE  BOTH  LARGE 


This  integral  appears  in  equation  (3.23).  We  shall  show  that  it  can  be  neglected 
when  compared  to  one  in  which  one  of  the  k’s  is  small.  For  this  purpose  it  is  desirable 
to  go  backwards  and  rewrite  the  required  integral  as  (we  shall  be  interested  only  in 
the  case  of  even  values  of  m) 


(1-2) 


71  +  72  —  71  COS  01  —  72  COS  02 


d0id02  • 


We  shall  not  be  too  rigorous  in  the  following,  but  this  as  well  as  several  integrations 
of  later  appendices  can  be  put  on  a  more  rigorous  basis  by  following  a  method 
developed  by  Duffin  [  17  ] .  When  the  m’s  are  large  the  integrand  in  the  neighborhood 
of  (0i,  02)  =  0  gives  the  main  contribution  to  the  integral.  Hence  if  we  introduce 
the  factor  J0(aR)  with  a  very  small  and  R  =  i<t>\y\/2  +  j<foy\/2  or  R2  =  027X  +  0272 , 
replace  the  denominator  of  the  integrand  by  |yi0i  +  £7202,  allow  the  integration 
to  extend  over  the  entire  real  (0 1,  0 2)  space,  and  transform  to  polar  coordinates, 
our  integral  becomes 


G-3) 


1  f 00  2RJ0(aR) 

2tt2(7i72)1/2  Jo  R 2 


e 


\Rk  cos  0 


dd  — 


iRk  '  cob  0  ' 

e 


where  k  and  k'  are  the  vectors 


k  =  wiry! 1/2  %  +  men1'2  j  ,  kr  =  m3 71 1/2  i  -p  me/21/2  j  , 

and  0  and  0'  are  the  polar  angles  between  the  vectors  R  and  k  and  R  and  k'.  The 
integrals  with  respect  to  0  are  again  Bessel  functions  of  order  zero.  The  desired 
expression  becomes 

(W)  f  -  Uk’mR 


-  if  ^  »  -  j^dR -r^[i- jmR }  • 

Since  both  of  the  integrals  in  the  bracket  are  “Bateman  integrals”  (see  Watson 
[22]),  the  entire  expression  becomes 

d-5)  =  P°8  *='/«  -  log  */«]  =  T{y^yn  lo6  <*'/*) 


1 

x(7i72)1/2 


log 


72m2  +  72mn 
SY2m\  +  yjmy  * 
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a  result  independent  of  the  value  of  the  small  number  a  which  was  introduced  in  the 
integrating  factor  J0(aR). 

The  case  of  interest  for  equation  (3.24)  in  the  text  is  mi  =  m2  =  m3  =  2m  and 
m4  =  0.  Our  required  integral  reduces  to 

(I_6)  x(7i72)1/2  l0g  (72  +  73)  ’ 

which  for  fixed  71  and  72  is  small  compared  to  [x(7i72)1/2]-1  log  2m  as  m  ->  °o . 


APPENDIX  II.  GENERALIZATION  OF  AN  INTEGRAL  OF  WATSON 


Let 

(Il-la) 

CO 

1 

II 

(Il-lb) 

-i/ir 

0 


d0id020d3 

(2  +  a2)  —  cos  —  cos  02  —  a2  cos  03 


d0id02d03 

sin2  |0 1  +  sin2  \<t>i  +  a2  sin2  ^03  ' 


The  special  case  a  =  1  has  been  discussed  by  Watson  [21].  We  shall  follow  his 
technique  here.  If  one  introduces  a  set  of  Cartesian  variables  Xi  =  tan  |0i,  x2  = 
tan  502,  and  x3  —  a  tan  §03,  then  transforms  these  to  spherical  polar  coordinates, 
and  finally  replaces  the  angle  0  by  0  =  5^,  the  integral  (Il-lb)  reduces  to 


(H-2)  1(a) 

_4 a  Pr  PT  _ sin  6  drdddft _ 

7rV0  Jo  Jo  a2+r2sin20[5a2sin20sin2'4r+(l+a2)cos20]  +  lr4(2+a2)sin40cos20sin2^r‘ 


The  integration  with  respect  to  6  can  be  carried  out  in  an  elementary  manner  if 
r  is  replaced  by  a  new  variable  t  which  is  defined  as  £21/2  =  r  sin  6  (for  fixed  0)  and 
the  order  of  integration  is  interchanged.  Then 


(II-3)  7(a) 

in* 


dd 


0  Vo  J o  a2  +  £2[a2sin2 0sin2^ +2(1  +a2)  cos2 0]  +  (2 +a2) £ 4cos2 fein2^ 

2  \/2  f~  P/2  dtd* 


ff 


V(l  +  <2  sin2  ¥)[a2  +  2t2(l  +  a2)  +  (2+a2)*4  sin  2*] 


The  integral  over  t  can  be  reduced  to  an  elliptic  integral  if  we  replace  ^  by  a  new 
variable  f  which  we  define  by  tan1!'  =  f/(l  +  t2)112. 

Again  we  interchange  the  order  of  integration.  We  find 
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(I I“4)  1(a) 

=  2s/2  7T-2  r 

Jo 

2a'“ 


dr 


f. 


dt 


a +r2)1/2 Jo  {(i+^)[a2(i+r2)+«2(r2(2+a2)+2(i+a2))]}1/2 


23/2  r° 
lr*  Jo 


dt 


’o  (l  +  r)1/2[r2(2+a2)+2(l+a2)]1/2 
The  transformation  f  =  tan  x  leads  to 

i.).er  — *■ 

a7T2  Jo 


(-[ 


1+r2 


r2(2+a2)+2(l  + 


“  l/2\ 

^)J  /  • 


(II-5)  7(a)  = 


K' 


1 


[20±^)_sin2x]1/2  ^[^)-sin2x]1/2; 


If  we  substitute 

a**)  *™--\£Z^rVA¥S'L 

into  (II-5)  and  interchange  the  order  of  integration  and  summation,  we  find 


(II-7)  1(a)  =  - 


‘HI  \d  f,  r»(n+4+«)  r° 
1 nr3  \_dt  n„o  w  !r(?i-j- 1  -(-2c)  J 0 


dx 


[2(1  +a2)a-2  —sin2  x]n+1  /2+«  J  e„0 


We  perform  this  last  integration  by  choosing  c  to  be  the  smallest  root  of 


(1 1-8)  c2  -  27c  -f  1  =  0 

with 

(1 1-9)  7  =  (4  +  3a2)/ a2  . 

Then 

(11-10)  4c |2^  a  ^  —  sin2  x|  =  (1  +  cc2lx)(l  +  ce_2,x)  . 


If  we  substitute  this  expression  into  (II-7)  and  expand  the  integral  in  a  power 
series  in  c  we  find 


(11-11) 


JV/2 
0 


[2(1  +  a)a  2  —  sin2  x] 


— n— X/2—  e 


dx 


=  |tt(4c) 


n+  1/2  +  e 


2F1  (n  -j-  \  +  e,  n  +  ^  +  e;  1 ;  c2) 


where  27\  is  the  hypergeometric  function  of  the  four  arguments  given  in  the  paren¬ 
theses. 

The  integral  (II-7)  can  be  expressed  in  terms  of  a  hypergeometric  function  of  two 
variables  if  (11-11)  is  substituted  into  (II-7)  and  the  definition  of  $4  is  applied.  We 
obtain 


(II-12)  1(a)  = 


l_/d^[(4c)1/2+«  r2(l  +  *) 


2l/2a7r2  l de 


r(i +2c) 


$i(h~\-e,  i+e;  l+2e,  1;  4c, c2)] 
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where  denotes  the  fourth  kind  of  Appell’s  hypergeometrie  function  of  two  vari¬ 
ables  [  see  equation  (4.20)  ] .  It  has  been  shown  by  Bailey  that  5U  can  be  factored 
into  a  product  of  two  ordinary  hypergeometrie  functions  of  one  variable  (see  also 
p.  238, [15]) 

(11-13)  5i[u,  v  +  v'  -  u  -  l;v,v';x(l  -  y),y(  1  -  x)] 


=  2Fi(u,  H -v'  —  u  —  1 ;  t>;  x)  2Fi  (u,  v  +  v'  —  u  —  1 ;v';y) 


provided  that  x  +  y  <  1  and 

(11-14)  {\x(l  -  y)|}v*  +  { |  7/(1  -  x)|}"2  <  1  . 

The  values  of  the  significant  parameters  in  our  case  are 

(11-15)  u  =  i  +  e 

(11-16)  v  =  1  +  2«,  v’  =  1 

(11-17)  v  +  v'  -u-l  =  \  +  e  =  a 

(11-18)  x(l  —  y)  =  4c  while  7/(1  —  x)  =  c2. 

The  explicit  values  of  x  and  y  are  the  roots  of  this  pair  of  equations  for  which 
(11-14)  is  satisfied.  The  value  of  c,  the  positive  root  of  (II-8),  is 

(11-19)  c  —  y  —  (y2  —  1)1/2 . 

Since  the  square  root  of  c  is 

(11-20)  c1'2  =  2-u*{(y  +  l)1'2  -  (7  -  1)1/2} 

we  see  that 

(11-21)  {|x(l  -7/)!}1/2+  {\yd  -Z)!}1'2-  1  =2c1'2  +  c-  1 

=  21'2{(t  +  1)1/2  -  (T  -  1)1/2}  +  (7  -  l)1/2{ (7  -  1)1/2  -  (7  +  1)1/2} 

=  -  [  (7  -  1)1/2  -  21/21  [  (7  +  l)1'2  -  (7  -  1)1/21 . 

As  long  as  a  is  real  [see  equation  (II-9)  ]  7  =  3  +  4a-2  >  3,  and  (11-21)  is  negative 
as  is  required  by  (11-14). 

The  roots  of  (11-18)  with  the  property  x  +  y  <  1  are  x  =  1  —  k\  =  ( kr3 )2  and 
y  —  k\  =  1  —  (fcg)2,  where 

(II-22)  h  =  my  ~  D1/2  +  (7  -  3)1/2]  [(7  +  1  -  (7  -  1)1/2] 

(11-23)  fe  =  §[(7  ~  1)1/2  -  (7  -  3)1/2]  [  (7  +  l)1'2  -  (7  -  1)1/2]  . 

Our  required  integral  then  takes  the  form 


(11-24)  1(a) 


1 


21/2cor2 


•  ||  [(4c)1/2+<  F( i+€,  i+e,  1+2.,  *;>»+«,  i+«,  1,  *•’)]} __0  ■ 

The  quantity  inside  the  bracket  is  of  the  same  general  form  as  the  corresponding 
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quantity  treated  by  Watson  in  his  analysis  of  (II-l)  with  a  =  1.  We  employ  his 
result  and  find 


(11-25)  1(a)  =  23/2k2k3K(k2)K(k3)/aT2 

where  K(k)  is  the  complete  elliptic  integral  of  the  second  kind.  Various  relations 
are  easily  shown  to  exist  between  k2,  ks,  k'2,  and  k[.  For  example, 


(11-26) 


k2k3  =  2  (k2k3)x/2 


(k3/k2)  =  \(y/y  —  3  +  y/y  —  l)2  . 


Equation  (11-13)  reduces  to  Watson’s  result  when  a  =  1  (and  hence  y  =  7)  for 
then  k2  and  k3  agree  with  Watson's  value  for  these  quantities. 

Equation  (11-25)  can  also  be  written  as 


(11-27)  1(a)  =  4[  (7  +  1)1/2  -  (y  -  l)ll2]7r~2a~1K(k2)K(k3)  . 


Two  limiting  forms  exist  for  I  (a)  in  the  range  of  very  small  and  large  a.  As  a  -*■  0, 


(11-28)  1(a)  ~  -  log  (26/2  a'1)  , 

TT 


a 

Figure  13 

Variation  of  the  generalized  Watson  integral  (equation  Il-la)  with  o. 
The  Watson  value  0.505  corresponds  to  a  =  1. 


while  as  a  — ►  <» 

o5/2/<-jl/2  _ 

(11-29)  1(a)  ~  - -1  [K( 21/2  -  l)}2  =  0.633a"1  . 

tt  a 

We  have  plotted  1(a)  as  a  function  of  a  in  figure  13. 
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The  function  1(a)  is  related  to  the  integral  (4.28)  through 


-iff! 


<201^02^03  -1  r/r  /  n  l/2\ 

- - - - - - —  =  72  /(LY1/Y2J  ) 

(272+71)— 71  COS  01— 72  COS  02— 72  COS  <^>3 


APPENDIX  III.  ASYMPTOTIC  FORM  OF 


(III-l) 


biis-- 


(m101+m^2+ro^3)  d0id02d0 3 


Z  7/(1  -  COS  0y) 

l 

AS 


2  2  —1  2  —1  2  _ 1 

s  =  mi7i  +  WI272  +  m3y3  — »  00 

This  type  of  integral  has  been  evaluated  by  Duffin  [17]  when  71  =  72  =  73.  We 
note  that  when  k  is  large  the  main  contribution  comes  in  the  range  of  small  values 
of  the  0’s.  The  integration  can,  without  significant  change,  be  extended  over  the 
entire  range  of  positive  0’s.  One  introduces  a  radius  vector  R  —  i<J>iy\n  +  j<hy\12  + 
&037*/2  so  that  the  polar  coordinate  representation  of  the  integral  becomes 

(III-2)  I  •  f  f  ~  ‘sin  mde  • 

where  s  is  the  vector  s  =  imiy~l/2  +  jm 272 1/2  +  km3y~l/2,  and  9  is  the  polar  angle 
between  R  and  s.  After  integration  with  respect  to  9  we  have 


(1 1 1-3) 
with 


1  f-  sin  Rs  _ 1 _ 

tt2(7i7273)1/2  Jo  Rs  2^(717273) 1/2 


1 1 1-4) 


/  2  -1  .  2  -1  -  2  — l\l/2 

s  =  (mi7i  +  m2 72  +  m3y3  ) 


This  asymptotic  form  to  (III-l)  is  equivalent  to 

(III-5)  e  (71  72  73  Iml(xyi)Im2(xy2)Im3(xy3)dx  ~  2irs(y  17273) 1/2 

as  s  -►  a> .  Here  the  Z’s  are  Bessel  functions  of  purely  imaginary  argument. 


APPENDIX  IV.  ON  THE  INTEGRATION  OF 


x-l/2  e-*(rl+y2)  Ia(xyi)Ip(xy2)dx 

AS 


(a2y  1 1  +  /S272  *)  — >  w 


IV-1) 


Fa.fi( 71,  72) 
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We  use  the  double  integral  representation  of  F 


(IV-2) 


1  rr _ 6’Wi^)  cut>  _ 

4ir 3/2  JJ  [71(1  —  cos  <£i)  -f-  72(1  —  cos<£2)]1/2  ' 

—  TC 


When  a  and  0  are  large  the  main  contribution  to  the  integral  comes  from  the  region 
of  small  $  1  and  <fa.  Hence  if  we  introduce  the  integrating  factor  exp  ( —  n R)  (here  n  is 
small  number),  with  R 2  =  7i<£2  +  7202,  we  can  extend  the  range  of  integration  over 
the  entire  real  plane.  We  then  have,  after  a  reduction  to  polar  coordinates, 

(IV-4)  Fa  A 71,  72)  ~  I  I"  ^_1{exp  (-/iB  +  cos  0)}  2£dfld0  , 

where  A:2  =  a27^  + 

Hence 

(IV-4)  ~  5  (— )’/!  f”  e-'V„(fcB)d«  =  (2,rv.72)“1,,0.!+^)',,!  • 

z  \'ryi72/  »/o 

In  the  limit  of  large  /c  and  small  n  we  find 

(IV-5)  F.,*(7fc  72)  -  (2ir)-1/2(o!272  +  /327i)-1/2  . 

This  result  could  also  have  been  obtained  by  applying  various  asymptotic  forms 
for  Legendre  functions  to  (4.41). 

o  -o  -o  -o-  o- 

In  conclusion,  the  author  wishes  to  thank  Mr.  J.  Bradley  for  his  aid  with  the 
calculations  that  were  required  for  the  construction  of  the  various  figures. 
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